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1. Introduction 

It has been well established that certain kinds of recursive T matrix algorithms (known 
as RCTMA) ^'^ are numerically stable and can be used to solve the Foldy-Lax multiple- 
scattering equations for particles exhibiting "modest" inter-particle couplings. By "modest 
couplings" , we refer to situations in which the order of orbital number of the Vector Spherical 
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Wave Functions (VSWFs) necessary to describe the field scattered by each particle in an 
aggregate of particles are not too much larger than that necessary for describing isolated 
particles. The "modest coupling" criteria apply to a host of multiple scattering situations, 
including systems of dielectric particles comparable in size to the wavelength and for most 
packing fractions including dense packing. The modest coupling criteria can also apply to 
metallic particles under certain conditions. 

Like any multiple scattering technique not employing matrix balancing, the RCTMA can 
encounter numerical difficulties in certain extreme situations of strongly coupled resonant 
phenomenon. In this work, we present a matrix balanced form of the Recursive Centered T 
Matrix Algorithm (or RCTMA) that can readily be employed even in the presence of strong 
{i.e. resonant) inter-particle couplings. The rather extreme situation of "strong couplings" 
studied here will generally require carefully micro-scaled engineered systems where high 
Q-factor resonances can occur for particles illuminated in isolation, and in which the parti- 
cles are sufficiently closely spaced that neighboring particles modify the resonance response 
properties. Examples of strong inter-particle couplings can be found in particles exhibiting 
plasmon resonances, surface resonances, or even photonic jet phenomenon. 

In section 2, the notation is introduced in a brief review of the relevant multiple-scattering 
theory. Section 3 describes an analytic matrix balancing procedure used to 'well-condition' 
the multiple scattering system of equations. A matrix balanced RCTMA is derived in section 
4. Essential formulas for applications are summarized in section 5. Their applications are 
then demonstrated by applying matrix balanced RCTMA calculations to study systems of 
interacting localized plasmon excitations. Some known and novel aspects of interacting 
localized plasmon excitations are presented herein. 

2. Multiple-scattering theory - VSWF approach 

Let us consider an arbitrary incident electromagnetic field incident on a collection of three- 
dimensional particles (as shown in fig.l). The particles are considered as 'individual' scat- 
terers if they can be placed in a circumscribing sphere lying entirely within the homogeneous 
medium (actually this constraint can frequently be relaxed, cf."). 

The electromagnetic field incident on an A^-particle system, Ei, is developed in terms 
of the transverse regular VSWFs developed about some point O arbitrarily chosen as the 
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Fig. 1. Schematic of an field incident on a collection of scatterers centered on xi,X2, ...,XAr. The radii of 
the respective circumscribing spheres are denoted Ri, R2, i?jv- 

system 'origin': 

00 n 

Ei (r) = ^0 X] X] {Rg [M„m (kr)] ai,„_^ + Rg [N„„ (A;r)] a2,n,m} 

n=l m=—n 
2 00 

= ^oJ2Y.^9 (kr)] = Eo Rg (A;r)] a (1) 

q=l p=l 

where Eq is a real parameter determining the incident field amplitude. Eq.(l), introduces a 
condensed notation for the VSWFs, M„„i and N„„i: ^'i^p (A;r) = M„,m (A;r) and ^'2,^ (kr) = 
Nn'm{kr). The notation Rg[] stands for "the regular part of and distinguishes these 
regular VSWFs from the "irregular" scattered VSWFs (cf. appendix A). In the second line 
of eq.(l), the two subscripts (n, m) are replaced by a single subscript p defined such that 
p{n,m) = n{n + 1) — m and has the inverse relations': 

n{p) = Inty/p m{p) = —p + n{n + 1) . (2) 

The last line of eq.(l), adopts the compact matrix notation allowing the suppression of the 
summation symbols. ' The superscripted ( * ) stands for the transpose of a column 'matrix' 
of composed of VSWFs into a row 'matrix' of these functions. 

For points external to all individual circumscribing spheres, the total field, Et (r) can 
be written as the sum of the incident field, and a set of 'individual' scattered fields, Ei"'\ 
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centered respectively on each of the particle centers: 

(r) = Ei (r) + E(^-) (r,) 

N 

= Eo Rg (A;r)] a + EoJ^^' i^r^) f^^ (3) 

i=i 

where each scattered field, Eg"''', is developed, with coefficients f^^\ on a basis of outgoing 
VSWFs defined with respect to the associated particle center, denoted Xj. The spherical 
coordinates relative to each scatterer are denoted r j = r — Xj . 

The crucial idea of Foldy-Lax multiple-scattering theory is that there exists an excitation 
field, Ecil (r) , associated with each particle which is the superposition of the incident field 
and the field scattered by all the other particles in the system (excluding the field scattered 
by the particle itself).'' From this definition, the excitation field of the j*'* particle can be 
written 

N 

E(£ (r,) = Eo Rg (A:r,)] ^ = E^ (r) + ^1'^ (r^ 

N 

J^^'^^a + J2 ^^'''^ fN 

where e^'' are the coefficients of the excitation field in a regular VSWF basis centered on the 
j*'^ particle. In the second line of eq.(4), we have used the translation-addition theorem^ ''^'■'^) 
and introduced the notation where J^^'^^ = J (kxj) is a regular translation matrix and 
, — Xi)] is an irregular translation matrix. Analytical expressions for the 
matrix elements of J {k:s.j) and H {k'^j) are given in refs.^'''. 

The other key idea of multiple scattering theory is that the field scattered by the object, 
is obtained from the excitation field e"^ via the 1-body T matrix, t[^\ derived when 
one considers the particle to be immersed in an infinite homogeneous medium. This relation 
is then expressed as 

/^■^ = T[^^e^!^ (5) 

[The index 1 on the t[^'^ indicates that this T matrix concerns an isolated particle, henceforth 
referred to as a '1-body' T matrix.] Employing eq.(5) in eq.(4), one obtains a Foldy-Lax set 



EoRg[^\kv, 



(4) 



5 



of equations for the excitation field coefficients^''': 

TV 

e!^) = J(^'°)a + H^''^ ^i'^ 4^ 3 = 1, N (6) 

For numerical applications where one is obliged to solve the equations on a truncated 
VSWF basis, it is advantageous to work with a set of formally equivalent equations involving 
the scattering coefficients This set of equations is derived by multiplying each of eqs.(6) 
from the left by Ti''^ and using eq.(5) to obtain 

TV 

In the RCTMA, one calculates the centered multiple scattering transition matrices, T^'^\ 
which directly yield the scattered field coefficients in terms of the field incident on the system 
through the expression 

N 

ff = Ti/''^) a^'^ with a(^-) ^ J(^'°) a (8) 

k=l 

In this equation, we have introduced the column matrix a^^'^ which contains the coefficients 
of the field incident on the entire system developed on a VSWF basis centered on the j*'* 
particle. 

3. Basis set truncation and matrix balancing 

Although the multiple scattering formulas of the previous section are expressed as matrix 
equations on VSWF basis sets of infinite dimension, the finite size of the scatterers natu- 
rally restricts the dimension of the dominate VSWF contributions. In order to discuss this 
phenomenon analytically, we consider the case of spherical scatterers. For non-spherical 
scatterers, the matrix balancing procedure described below should be applied to the circum- 
scribing spheres of the particles. 

The Mie solution for a sphere of radius Rj immersed in a homogeneous host medium, can 
be cast in the form of a 1-body T matrix that is diagonal in a VSWF basis centered on the 
particle: 

= Sg,g'5pyTi{j,n{p),q) (9) 

<?,p;<?',p' 

where the T^'" {n{p), q) correspond to the Mie coefficients and depend on q and n (cf. eq.(l)). 
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With the objective of matrix balancing, it is helpful to express the Mie coefficients of the 
scatterers in terms of the Ricatti Bessel and Hankel functions, respectively ipn [z) = zjn {z) 
and (z) = zhn {z), and their logarithmic derivatives 

$ (z) = ^ (z) = do) 

" - {z) - a (z) ^'^^ 

The T matrix elements of eq.(9) for a sphere of dielectric contrast pj = kj/k can then be 
cast in the convenient form' : 

T(.. ^^ jkR,) 'f^n jkR,) - P,<^n {p,kR,) _ jkR,) 

' en (kR,) p,$„, (p.kRj) - ^vl>„ (kR,) - U (kR,) ' 

Tf7 n 2) = 'f^n {p,kR,)-p,<!>^{kR,) ^ i^^jkR,) 

^^'""'^^ en(A:i?,) p^^r,{kR,)-f<!>r,ip,kR,) - U{kR,) ^^''''^^ ^''^ 

where k is the wavenumber in the external medium. The normalized T matrix coefficients, 
T{j,n,q), of eq.(MieT) contain a rich resonant structure. The ratios ifjnikRj) /^nikRj) 
on the other hand have an exponentially decreasing behavior for large, n ^ kRj as is 
demonstrated in fig.2 for kR = 10. One can remark from figure 2 that \ipn (kR) /^n {kR)\ 
become quite small beyond nmax = kR + 3 and its value at n = 14 is ~ 2 10~^. Although 
these factors permit an appropriately truncated VSWF basis set to contain essentially all 
the physical information necessary for accurate calculations, they also tend to produce ill- 
conditioned linear systems when one is obliged to enlarge the VSWF space far beyond 
~ kR + 3 in order to account for strong coupling phenomenon. 

A solution to the above problem is to 'balance' the matrix manipulations in section 4 
below by defining "normalized" scattering and incident coefficients: 



n,,^^ni,)ikR,)[a(\^ (12) 

For notational purposes, it is convenient to define diagonal matrices [ip^-'^ and [e'--'^] with 
Ricatti-Bessel functions along their diagonals, namely [''P^''^]^, p,.qp = Sq,q'Sp,p'i'n{p) (kRj) and 
[^^''^Iq' p'-qp ^ ^q,q'^p,p'^n{p) (kRj). This uotatiou allows normalized or 'balanced' versions of 
the one-body and many-body T-matrices to be defined respectively as 

^0-) ^ ^^0-)] tU) [^0-)] -1 and T^^ ^ [e^^T T^^ [^^'^] (13) 
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Fig. 2. Plot of the spherical Bcsscl to Hankcl function ratio, ipn (kR) (kR) occurring in the Mie 
coefRcients when kR ~ 10. 

In terms of these normalized quantities, eq.(8) then reads 

N 

fS = Y.'^f'^''' J = l,...,iV (14) 

k=l 

In the next section, these 'normalized' T^"*^ and T^^^^"^ are used to derive a matrix balanced 
version of the recursive T matrix algorithm. 

4. Derivation of a matrix balanced recursive algorithm 

In this section, we derive a matrix balanced version of the Recursive Centered T Matrix 
Algorithm (RCTMA) using purely algebraic manipulations. The recursive algorithm can be 
invoked once we have a solution for the Tjy'_i matrices of a > 1 particle system. If we 
wish to solely use the recursive algorithm to solve a system, we initiate the recursive process 
with a single particle solution described by t[^'^^ = T^^'. 

One then considers an arbitrarily positioned particle being added to the system. The 
excitation field on a particle A^ added to the system can be expressed as the superposition 
of three fields. The first contribution is simply the field incident on the system, the second 
contribution results from the scattering of the incident field by the A^ — 1 cluster of particles 
onto the particle A^, and finally the third contribution comes from field scattered by the 
particle A^ onto the A^ — 1 cluster and which returns to the A^*^ particle as an excitation 
field. Invoking the translation-addition theorem and eq.(8), these three contributions can 
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be expressed in matrix form as^ 



N-1 



IN 



(15) 



j,k=l j,k=l 

Defining now the normalized irregular translation matrices and excitation coefficients 
respectively as 



H^^'^) ^ [^0-)] [^(fc)] -1 and e^^ ^ [^(^'T e^^} 



(16) 



the normalized form of eq.(15) can be written 



j,k=l j,k=l 



(17) 



where we also used the definitions in eqs.(12) and (13). 

Recalling that the excitation field is linked to the scattered field by the 1-body T-matrix 
via eq.(5), and invoking the definitions of eq.(12) and (16) we can write 



f 



N 



(18) 



Employing this relation for particle N on the LHS of eq.(17) and rearranging we obtain 

N-l 



—{N,j)^{j,k) ■jj(k,N) [ -j{N) 

In 



j,k=i 

N-l 

j,k=l 



(19) 



— (N N) 

We now take the normalized ' matrix to be expressed as 



N 



-1 



^ ■Jj{NJ)7j^{j,k) ■j^(k,N) 
j,k=l 

-(N,N) 



(20) 



With this assignment, we multiply both sides of eq.(19) by ' and obtain an expression 
consistent with equation (14): 

N-l 

j{N) _ -{N,N) (^) -{N,N) —{N,j) -{j,k) (fc) 

Jn —-'-N " + N 2-^^ -'-N-l"' 

3,1=1 



N-l N 



(21) 



1=1 



1=1 
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(^N,k) 

where we have assigned the matrix ' , k N a.s 

N-l 

^ N = J- N ^ N-l l^^j 

i=i 

One completes the description of the scattering by the system by remarking that the field 
scattered by the other particles in the system are the superposition of the field that would 
be scattered by the — 1 cluster in the absence of the A^*^ particle, plus the field scattered 
from the — 1 particle originating as a scattered field emanating from the A^*^ particle. 
Using again the translation-addition theorem, the field coefficients of can in turn be 
expressed in a form consistent with equation (14) as 

N-l N-l 



In = 2^ ^ N-iO' + 2^ ^ N-1^ In 

l=k k=l 

N-l N-l 

ETp{j,k) (fc) \^ 7p{j,k) -Yrik,N) 7p{N,N) (jy) 

1=1 k=l 
N-l N-l 

+ Y.T.'^NW^'''T'N^''ai'' 

1=1 k=l 

N-l N 

= t5^'^^ a(^) + J2 Tn'^ = E ^^n'^ (23) 



k=l k=l 



where we invoked eq.(21). In the second and third lines we have defined the T^' and T^' 



matrices such that 



Af-l 

Tn =2^ ^iv-i H (24a) 



k=l 

N-l 



7r(3,k) 7j^(hk) , sr^Tj^[],l) -jj{l,N)j^{N,k] , 
J- N = J TV-l + 2^ J 7V-1 ^ J- N l^4bj 

1=1 

At this point, all the T ^ matrices have been obtained and the matrix manipulations 
in eqs.(20), (22) and (24) can then be repeated to add as many particles to the system as 
desired. 

A. Relationship with system matrix inversions 

Although the recursive algorithm is quite efficient for systems with relatively small numbers 
of particles, for systems with many particles, one may prefer to try and solve an entire 
A^-particle system directly. A balanced linear system for the entire system corresponding to 
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our recursive algorithm can be obtained by applying the relation of eq.(5) to the left hand 
side of eq.(7), then multiplying both sides of the resulting equations by the [ip^-'^ matrix 
and finally rearranging to obtain a system of balanced linear equations for the unknown 
scattering coefficients: 

N 



f 



\^ T7(J»-f(*^) -(7) 



N 



J 



(25) 



where we used eqs.(12) and (13). The system of linear equations in eq.(25) can in principle 
be directly solved by inverting the balanced system matrix: 

-1 



' 7(1) ' 

In 
7(2) 
In 




In 





(1) 



-H 



(1,2) 



-H 



{1,7V) 



-H 



<2,1) 



(2)- 



-H 



{2,7V) 



-1 

















(26) 



Once we have inverted this system, one can associate each block with a corresponding 
matrix as shown below as 



' 7(1) " 

In 




7(2) 

J N 




-riN) 

In 





7p(l,l) 7p(l,2) 
-'7V Af 

7^(2,1) ;=(2,2) 
-'7V N 



^{7V,1) -(7V,2) 
-'7V N 



7p(l,^) 
N 

(2,7V) 



T 



N 



7f^(N,N) 
N 



7(1) 



(27) 



which is the same form as the desired solutions given in eq.(14). 

5. Summary and applications to localized plasmon excitations 

In this section, we will apply the RCTMA to solve systems exhibiting strong interactions 
between localized plasma resonances. We begin this section by summarizing the balanced 
recursive algorithm. We then recall some useful formulas for extracting physical quanti- 
ties from the T matrix. Finally, we carry out some illustrative calculations for strongly 
interacting systems. 



A. Summary of the balanced RCTMA algorithm 



(1) T.(2) 



In order to implement the RCTMA, one must first solve the 1-body T- matrices, T{ ,T{ 



(7Vtot) 



for all the particles in the system. Normalized versions of the 1-body T matrices 



11 



and the irregular translation matrices^'', H^^'^\ are then calculated via 

rf ^ [e^^] t[^^ [^(^-^j F^^''^ = [^(^-^j H^^'^^ [^e^)] (28) 

where the diagonal matrices [V^^^^] ^ = 5q,q>5p,p4n{p){,kRj) and [^^^^] ^ ^ p, = 

^q,q'^p,p'^n{p) i^Rj) respectively have Ricatti-Bessel and Ricatti-Hankel functions on the di- 
agonal. {Rj being the radius of the circumscribing sphere of the jth scatterer). 

The balanced recursive algorithm is that the solution for the T^'^^ matrix is obtained 
from the T matrices of the — 1 system, T^'^l, via the matrix inversion in eq.(20). All 
the other matrices T^' with j ^ N or k ^ N are then obtained via matrix multiplications 
and additions via equations (22) and (24). This process is then repeated as many times as 
desired. 

B. Physical quantities 

When the incident field is a plane wave, it is convenient to express physical quantities in 
terms of cross sections. Appealing to the far-field approximation of the field, the extinction 
and scattering cross sections of clusters of A^ objects can be respectively expressed*^'^ 



c^ext = --^ Re 



N 



and a,.., = ^ /^'^'V^^'^Vi'^ (29) 



fc2 

j,fc=l 



.k=l 

It is also possible to produce analytical expression for local field quantities like indi- 
vidual absorption cross sections. For lossy scatterers in a lossless host medium, one can 
obtain individual particle absorption cross sections by integrating the Poynting vector on a 
circumscribing sphere surrounding the particle to obtain the formula as 



cjij) - -_Rei f(^)'teO) , 



(i) 

N 



(30) 



In an analogous fashion, optical forces on the particles can be calculated by integrating 
the Maxwell tensor on a circumscribing sphere surrounding the particle^"'^^. It is frequently 
convenient to characterize the optical force by vector 'cross sections', (Topt; defined such 
that the time averaged optical force on particles immersed in a liquid dielectric of refraction 
index nmed can be expressed as 

Fopt 1 1 Sine 1 1 ^ opt 

(31) 
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where ||Sinc|| = H^R-elE*^^ x Hinc}|| is the incident irradiance. The binding force and its 
associated cross section, o"b; between two particles separated by a relative position vector 
Tpos = r2 — ri, can be defined as 

F^^\ (F2 - Fi) -^pos = llSincll (32) 

C. Interacting localized plasmon excitations 

For those conductors, such as the noble metals, that support surface plasmon resonances, 
one can usually observe localized plasmon resonances in sufficiently small particles. These 
resonances are typically dominated by absorption if the particles are sufficiently small with 
respect to the incident wavelength and by scattering for larger particles. We chose to study 
silver spheres 50 nm in diameter immersed in air (for which both scattering and absorption 
are non-negligible). 

The study is carried out for wavelengths ranging from the near ultra-violet through the 
visible (300 to 850 nm). We ignore the relatively modest finite size corrections to damping^^ 
and simply adopt the bulk dielectric constant of silver from ref.^'^ and extrapolate between 
the experimental values. The extinction, scattering and absorption cross sections for these 
particles are readily obtained from Mie theory and are displayed in figure 3 as a function of 
frequency. These spheres are quite small with respect to visible wavelengths, (size parameters 
in the 300<-i>800 nm wavelength range go through kR = 0.52 0.20) and the isolated particle 
cross sections are obtained to high precision with rimax = 4. One can also see from figure 3 
that the strength of the plasmon resonance for these particles is about half due to absorption 
and about half due to scattering. 

One of the principal sources of interest of the localized plasmon resonances is their ca- 
pacity to produce large field enhancements in regions much smaller than the incident field 
wavelength. This property is demonstrated in fig. 4a) with a 2D and ID plot of the electric 
field intensity in and near an isolated 50 nm diameter silver sphere illuminated near its 
resonance peak (Aq = 365 nm with A^Ag = 0.077 + 1.6i). The plots in Fig.4 are performed 
in a plane containing the center of the sphere and perpendicular to kinc (the polarization 
lies along the horizontal axis). The dimensionless extinction and scattering 'efficiencies', 
Q = a/{nR^), at this frequency are are respectively Qcxt = 14.48 and Qscat = 6.76. 

We now use the balanced recursive technique to calculate the optical response of a dimer 
composed of 50 nm diameter silver spheres whose surfaces are separated by 1 nm. Although 
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Fig. 3. Total cross section 'efficiencies', Q = (j/{tiR^) for an isolated 50 nm diameter sphere. 

the T matrix calculated by RCTMA contains information for arbitrary incident fields, we 
study the physically interesting case of a plane wave perpendicular to the axis separating 
the particles. As is widely known, the response then depends strongly on the polarization 
of the incident light. In figures 5a) and 5c), the extinction and scattering cross sections 
per particle are presented when the polarization is respectively perpendicular and parallel 
to the symmetry axis. From figure 5c), one sees that the cross section for the parallel to 
axis polarization presents a two sphere coupled resonance that is strongly red-shifted with 
respect to the isolated particle resonance. The polarization perpendicular to this axis on 
the other hand presents only slight modifications with respect to an isolated sphere. 

The optical binding force cross sections for these same polarizations are respectively plot- 
ted in figures 5b) and 5d). While the binding force for the polarization perpendicular to the 
particle axis (cf.5c)) is slightly repulsive, the force for polarization parallel to the resonance 
can be highly attractive with the dimensionless \Qh\ attaining amplitudes of three orders 
of magnitude. There has already been experimental and theoretical evidence supporting 
the existence of optical force couplings in particles with plasmon excitations^ ^ although such 
high precision calculations at such small separations seems not to have been presented before 
now. 

This dimer system dramatically illustrates the 'strong' coupling category since correct 
calculations require that the VSWF space be enlarged far beyond the predominantly dipolar 
response characterizing the particles in isolation. The normalized cross sections per particle 
are given in the table 1 for different values of the VSWF space truncation. Although it was 
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Fig. 4. Electric field intensity ||Et|p/||Einc|P in an isolated 50 nm diameter sphere (Aq = 365 nm, 
^Ag = 0.077+ l.Qi). In a) is presented a 2D (hot) plot of the electric field intensity in a plane perpendicular 
to the wavevector and containing the origin of the sphere (the horizontal axis lies along the polarization 
direction). Fig b) is a ID plot of the field intensity along the line in this plane containing the direction of 
electric field polarization. 

necessary to go to ~ rimax = 30 to achieve 4 digit precision in all the cross sections, the table 
indicates that results were already quite good at rimax = 20. 

A base 10 logarithmic intensity field map for a two silver sphere dimer illuminated with 
light polarized along the symmetry axis (frequency near the coupled sphere resonance maxi- 
mum (Ao = 467 nm and A^Ag = 0.048 + 2.827i) is presented in figs. 6a) and figs. 6b) which are 
respectively a 2D plot (in the same plane as figure 4) and a ID plot along the symmetry axis. 
The size parameter of the individual spheres is kR = 0.34 and the isolated cross sections at 
this frequency are Qcxt = 0.136 and Qscat = 0.0963. As can be seen in fig. 6, the fact that one 
had to go so far beyond the dipolar response has a dramatic effect on the field inside and 
near the the particles. Notably, the fields inside the particles are no longer quasi- const ant 
as was the case for isolated particles. 
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Fig. 5. Dimensionless cross section 'efficiencies' per particle, Q = a/{2'KR^) and binding force 'efficiencies' 
for a dimcr of 50 nm diameter spheres (1 nm separation). In a) and b) the polarization is perpendicular to 
the symmetry axis and in c) and d) it is parallel to the symmetry axis. 



'^max 


5 


10 


15 


20 


25 


30 


35 


40 


Qcxt/2 


4.60 


15.53 


17.38 


17.20 


17.14 


17.13 


17.13 


17.13 


Qscat/2 


3.51 


10.62 


11.30 


11.04 


10.98 


10.97 


10.97 


10.97 


Qb 


-417 


-3639 


-5530 


-5918 


-6000 


-6015 


-6018 


-6018 



Table 1. Dimensionless cross section 'efficiencies' per particle in function of the VSWF truncation, ?imax- 
Qcxt = foxt/CSTri?^), Qscat = crscat/(27r_R^) and Qb = Cb/(7ri?^). The system is a dimer composed of 
D = 50 nm diameter silver spheres (1 nm separation) at (Aq ~ 467 nm and A^Ag ~ 0.048 + 2.827?) 

An important word of caution should be made at this point. Although 1 nm separation 
may appear to 'nearly' touching, the coupled resonance is in fact quite sensitive to exact 
separation details when resonant particles are so closely separated. For example, at a sep- 
aration distance of 0.5 nm for the silver dimer, the coupled plasmon resonance is displaced 
to Aq — 516 nm as compared with Aq — 467 nm for a 1 nm separation, and the multipole 
order has to be pushed to nmax — 50 to achieve four digit accuracy in the cross sections. 
Nanometer scale separations are not necessarily theoretical idealizations however as recent 
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Fig. 6. Logarithmic scale plots of the field intensity for a two sphere dimer with (Aq ~ 467 nm and incident 
light polarized along the sphere axis (A^Ag = 0.048 + 2.827i). In a) is a 2D plot in the plane containing the 
centers of the spheres and the polarization vector while b) is a ID Logarithmic plot along the symmetry 
axis of the spheres. Figures, c) and d) are the same as a) and b) respectively but for a 5 sphere chain of 
spheres at its resonance maximum (Ag = 561 nm and A^Ag = 0.0564 + 3.685i). (cf. figure 7). 
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experiments with DNA separators have demonstrated^'"'. Nevertheless, in apphcations hke 
DNA separators, one may well have to consider the strong optical forces between these 
particles on account of the exceptionally strong attractive optical forces efficiencies of these 
resonances. For instance, the binding force efficiency at 0.5 

mathrmnm separation was calculated at Qb = —20174 for Aq — 515.6 nm (cf. Qb = —6018 
for 1 nm separation at Aq — 467 nm). The question of perfect spheres exactly in contact 
however seems untenable from an experimental standpoint and quite difficult from a the- 
oretical standpoint on account of the singular behavior of the contact point. Theoretical 
'separations' of lA for instance require multipole truncations of the order of nmax ^ 120 
before convergence is achieved, but the idea of 'perfect' spheres separated by atomic scales 
has clearly gone beyond domain of applicability of our mesoscopique physical model in any 
case. 

It is also important to verify that the recursive algorithm works for more complicated 
systems. Towards this end, we illustrate in figure 7, the results of calculations for a system 
composed a line of 5 identical silver spheres spheres separated by 1 nm. For the binding 
force, we now present Qh,i which is the binding optical force between outermost spheres and 
their nearest neighbor and (5b,2 which is the binding force between central sphere and each of 
its nearest neighbors. It is interesting to remark that addition of other spheres in the chain 
dramatically lessens the strong binding force interactions between spheres even though the 
fields between the spheres (cf. fig.5) can still be almost as high as the dimer case. 

We remark that the interactions have continued to red-shift and widen the coupled 
"chain" resonance. This chain resonance peaks at ^ 561 nm and A^Ag — 0.0564 + 3.685i 
{Qext/5 = 14.416 and Qscat/5 = 12.543). It is clear from figure 7c) that the extinction 
cross section of the chain resonance is increasingly dominated by scattering rather than 
absorption. The number of VSWF orders necessary for high precision was also seen to de- 
crease slightly for the chain. The calculation of fig. 7 was carried out with Umax = 20 since 
calculations at rimax = 24 produced negligible differences on this scale. 

Despite the dominance of scattering, a considerable amount of absorption is still present 
in the 5 sphere chain. Furthermore, from the field maps in figures 6c) and 6d), one can see 
that the central sphere has the highest field internal field intensities, and one consequently 
expects increased absorption in the central sphere. This supposition can readily be confirmed 
quantitatively by using eq.(30) to calculate the absorption in each individual sphere. The 
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Fig. 7. Total cross section and binding 'efficiencies' for a line of 5 'touching' silver spheres 50 nm in diameter 
(1 nm separation). In a) and b) the polarization is perpendicular to the symmetry axis and in c) and d) it 
is parallel to the symmetry axis. 

results are given in table 2. 



Qa.,1 


Qa.,2 


Qa,3 


Qa.,4: 


Qa,5 


0.8346 


2.333 


3.030 


2.333 


0.8346 



Table 2. Individual absorption efficiencies Qa.j = faj/(7ri?^) in a five sphere chain at Ao = 561 nm and 
iVAg = 0.0564 + 3.685i 

We conclude this section with some calculations systems concerning larger chains of 
particles. One can remark that the chain coupled resonance continued to red-shift and 
widen when passing from the dimer to the five particle chain. Results for the extinction and 
scattering cross sections chains of 10 and 20 sphere chains are presented in figure 8 for the 
same polarizations and incident directions as considered previously. 

One readily sees that ultraviolet and perpendicular responses per particle seem to have 
stabilized for large chains. The collective chain response on the other hand continued to 
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broaden and slightly red shift as one passed from 10 to 20 sphere chains and it is an interesting 
point for future studies to examine the evolution of this phenomenon for even longer chains 
and to study the impact of defaults in the chains. 
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Fig. 8. Total extinction and scattering cross section efficiencies per particle in chains of 10 and 20 particles. 
In a) the polarization is perpendicular to the symmetry axis and in b) it is parallel to the symmetry axis. 

6. Conclusions 

The balanced recursive algorithm can give useful and highly accurate information in systems 
with large numbers of strongly interacting resonances. This has been demonstrated herein 
for the case of localized plasmon resonances and the studies presented here suggest that 
chains of closely spaced localized plasmons can have potentially interesting applications 
with respect to frequency shifting and broadening. Although not demonstrated here, this 
technique also proves useful for treating closely spaced systems possessing surface resonances 
of 'whispering gallery' type. 

It is worth remarking that matrix balancing seems to be a useful method to employ in 
almost any Foldy-Lax equation solution scheme, be that for direct system matrix inversion, 
iterative techniques or linear system solutions. In fact, some modern matrix inversion pro- 
grams actually integrate numerical matrix balancing into their algorithms. Nevertheless, 
since the matrix balancing in Foldy-Lax equations can be obtained analytically at relatively 
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low computational cost, it seems beneficial to carry out this balancing explicitly rather than 
relying on purely numerical treatments. 

The matrix balanced RCTMA has potentially interesting applications for other kinds of 
resonance phenomenon, notably whispering Gallery modes. Such studies are currently un- 
derway. Furthermore, the ability of the matrix balanced RCTMA to study defaults and small 
modifications in large complicated systems is particularly promising and will be employed 
in subsequent studies. 

Brian Stout and Alexis Devilez would like to thank Ross McPhedran, Evgeny Popov 
and Nicolas Bonod for helpful discussions. This work was funded in part by the grant 
ANR-07-PNANO-006-03 ANTARES of the French National Research Agency. 

Appendix A: Vector spherical wave functions 

The vector spherical wave functions can be readily written in terms of the Vector spherical 
harmonics (VSHs) and outgoing spherical Hankel functions : 

(kr) = Mr,Ukr) = ht (kr) X„^(^, 0) 
*2,p (kr) = N„^(A;r) = ^ [^n{n + l)ht (kr) Y^^e, 0) + [krh+ (kr)] ' Z„^(^, 0)] (Al) 

In the same manner, the regular VSWFs are obtained by replacing the spherical Hankel 
functions in eq.(Al) by spherical Bessel functions. Our adopted definition of the VSHs is 

Ynm{e, 0) = vYn^{e, 0) Z„„(^, 0) = . "'"^ X„„(e, 0) = Znm{e, 0) A? (A2) 

^^n[n + 1) 

where the Ynm{d^ 0) are the scalar spherical harmonics. 
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